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Abstract 

Thermodynamic property calculations of mixtures containing carbon dioxide (CO 2 ) and wa¬ 
ter, including brines, are essential in theoretical models of many natural and industrial processes. 
The properties of greatest practical interest are density, solubility, and enthalpy. Many models 
for density and solubility calculations have been presented in the literature, but there exists only 
one study, by Spycher and Pruess, that has compared theoretical molar enthalpy predictions 
with experimental data [T]. In this report, we recommend two different models for enthalpy cal¬ 
culations: the CPA equation of state by Li and Firoozabadi [2], and the CO 2 activity coefficient 
model by Duan and Sun [3]. We show that the CPA equation of state, which has been demon¬ 
strated to provide good agreement with density and solubility data, also accurately calculates 
molar enthalpies of pure CO 2 , pure water, and both C02-rich and aqueous (H20-rich) mixtures 
of the two species. It is applicable to a wider range of conditions than the Spycher and Pruess 
model. In aqueous sodium chloride (NaCl) mixtures, we show that Duan and Sun’s model yields 
accurate results for the partial molar enthalpy of CO 2 . It can be combined with another model 
for the brine enthalpy to calculate the molar enthalpy of H20-C02-NaCl mixtures. We conclude 
by explaining how the CPA equation of state may be modified to further improve agreement 
with experiments. This generalized CPA is the basis of our future work on this topic. 


1 Introduction 

Mixtures of carbon dioxide (CO 2 ) with water, including brines, are important in a number of 
industries, such as food products, combustion systems, cosmetics, and petrochemicals. They also 
play a major role in Earth science applications like geological carbon sequestration, in which CO 2 is 
injected and stored in saline aquifers located deep in the Earth’s subsurface la El El El El El [in]. In all 
of these applications, accurate models of thermodynamic properties are essential. The models must 
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often be reliable over a wide range of temperatures and pressures. This is the case, for example, 
in the work of the National Risk Assessment Partnership (NRAP), which is a collaborative effort 
between several United States Department of Energy national laboratories m- One of NRAP’s 
main missions is to predict the risk of CO 2 leakage from storage reservoirs (deep saline aquifers) to 
overlying shallow aquifers, which are potential sources of drinking water. The models employed by 
NRAP need to compute thermodynamic properties in the range of conditions between those of the 
shallow aquifers, which are roughly between 290-300 K and 3-20 bar, and those of the CO 2 storage 
reservoirs, which are between 323-423 K and 100-500 bar, depending on the subsurface depth. 

Thermodynamic properties of greatest practical interest can be classified into three categories: 
volumetric (e.g., density), phase behavior (e.g., solubility), and thermal (e.g., enthalpy). For mix¬ 
tures with CO 2 and water, there are several different models that can be used to compute the 
density and/or solubility over a wide range of conditions. Enthalpy models, however, are much 
more rare. This may be a reflection of the fact that enthalpy measurements are far less common 
than density or solubility measurements. To the best of our knowledge, there has been only one 
reference that has compared theoretical molar enthalpy (or specific enthalpy) predictions with ex¬ 
perimental data [T] . This study by Spycher and Pruess has been limited to pure components and to 
mixtures of CO 2 with pure water, not brines. Mixtures with pure water may be of direct interest to 
the industrial applications listed above. For carbon sequestration, fresh water is sometimes used as 
a surrogate for the brine in which CO 2 is dissolved. Nevertheless, in some cases it may be desirable 
to represent the saline environment of the aquifers with true brine properties. 

The main goal of this report is to identify models that can accurately calculate the enthalpy 
of mixtures containing CO 2 and water. For mixtures with pure water (not electrolytic solutions), 
we promote use of the cubic-plus-association (CPA) equation of state (EOS) developed by Li and 
Firoozabadi [2]. The CPA equation of state is applicable to pure components as well as to aque¬ 
ous mixtures containing substances such as CO 2 , hydrogen sulfide (H 2 S) and hydrocarbons. For 
these mixtures, it has been demonstrated to provide good agreement with experimental density 
and solubility data, but no comparisons have been made with enthalpy data. In fact, simultaneous 
prediction of phase equilibria and enthalpy is considered one of the remaining challenges for CPA 
equations of state m- In fluids where water is not present, such as in hydrocarbon-C02 mixtures, 
the CPA reduces to the widely-used Peng-Robinson equation of state. The CPA has been suc¬ 
cessfully implemented in codes that simulate three-phase (oil, water, gas) flows [TH]. For aqueous 
mixtures of CO 2 dissolved in brines, we recommend using the Duan and Sun activity coefficient 
model to compute the partial molar enthalpy of CO 2 [3]- We consider brines composed only of 
water and sodium chloride (NaCl), which is the predominant salt found in deep saline aquifers. 

The report is organized as follows. We present molar enthalpy calculations from the CPA 
equation of state for pure components and CO 2 -II 2 O mixtures in Sections and respectively. 
We show that it is applicable to a wider range of conditions than the Spycher and Pruess model [I] . 
In an effort to make the CPA EOS accessible and readily usable, extensive details- much more 
than presented in the original research paper [2]- are included in Appendices and Enthalpy 
calculations of H20-C02-NaCl mixtures using the Duan and Sun model are presented in Section [T2] 
and Appendix Although CPA accurately predicts molar enthalpies, it exhibits relatively poor 
agreement with experimental excess enthalpy data. Section together with |B.2f|B.4] explain how 
the Li and Firoozabadi version of the CPA equation of state may be generalized to improve the 
agreement with excess enthalpy data. This generalized model is the starting point of our future 
work that we outline in Section [5l 
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2 Enthalpy of pure components 


2.1 COo 


The molar enthalpy h is given by (301, which for a pure component fluid simplifies to 


h{T, P) = -RT^ 


5 in 0(r, P) 

df 




P,n 


( 1 ) 


The fugacity coefficient (j) of the pure component can be obtained from (851, which is the expression 


for the fugacity coefficient from Li and Firoozabadi’s version of the CPA EOS. Since CO 2 is not self¬ 
associating (hydrogen bonds are unable to form between CO 2 molecules), this expression simplifies 
tremendously because the association contribution to In cj), which involves the site fractions Xi: is 
zero. For these non-self-associating substances, Li and Firoozabadi’s CPA EOS is equivalent to the 
Peng-Robinson EOS. The molar enthalpy of CO 2 in the ideal gas state is found by integrating, 
with respect to temperature, the correlation for the constant pressure heat capacity presented in 
Appendix C of the textbook by Smith, Van Ness, and Abbott HU: 


/i« = iiLr + 6T_ A (2) 

where a = 5.457, b = 1.045 x 10“^ and c = —1.157 x 10® K^. For pure components, we 
compare theoretical calculations to experimental data from the National Institute of Standards 
and Technology (NIST) [12]. In order for the calculations to be compatible with the data, the 
same reference point for the enthalpy must be chosen. NIST presents their enthalpy data using a 
reference point where the internal energy of liquid water along the vapor-liquid (VLE) equilibrium 
curve at 273.16 K is assigned to be zero. With this choice of reference point, NIST reports the 
molar enthalpy of CO 2 at 273.16 K and 34.861 bar (the CO 2 vapor pressure at 273.16 K) to be 
8.804 kJ/mole. Therefore, the molar enthalpy at any other temperature T and pressure P is 


h(T, P) = -RT^ ( + psm - h(T = 273.16 K,P = 34.861 bar) + 8.804, (3) 

\ dT J p n 

where h{T = 273.16 K, P = 34.861 bar) is the right-hand side of 0 evaluated at the indicated 
conditions. Computing the enthalpy in this way circumvents the issue that the reference point for 
the CPA EOS, which is used to calculate the fugacity coefficient, may not be the same as that for 
the correlation in (j^. 

The CPA EOS (or more precisely, the Peng-Robinson EOS) produces accurate CO 2 enthalpy 
predictions over a wide range of temperature and pressure conditions (Figure [^, including those 
stated in Section for NRAP. The maximum error is less than 0.60 kJ/mole in virtually the 
entire range considered. Although not shown here, we have found that calculations for isotherms 
above 500 K follow the trend in Figure of increasingly better agreement with data at higher 
temperatures. The results are comparable to those from Spycher and Pruess, where the accuracy is 
stated to be within a few percent [T] . The parameters of the EOS have been fitted to experimental 
vapor pressure measurements, so that the vapor pressure predictions closely match the data, even 
near the critical temperature. Except for near the critical temperature, the CPA EOS can also 
accurately predict the molar enthalpy of vaporization A/ivap [Figure |^b)], which is defined as the 
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(b) 



Figure 1: CPA EOS-predicted enthalpies of pure CO 2 vs. pressure along five isotherms. Theoretical 
calculations are depicted in solid lines and are calculated from Q. Experimental data (dots) are 
obtained from NIST [12]. The absolute error of the predictions in (a) is shown in (b). 


difference in the vapor-phase and liquid-phase enthalpies shown in Eigure j^a). The enthalpy of 
vaporization rapidly approaches zero near the critical point; most of the error near this point comes 
from the liquid-phase enthalpy predictions. Due to thermodynamic singularities that arise, it is 
not uncommon for equations of state to have larger errors near the critical point. 


2.2 H 2 O 

The molar enthalpy of pure water is calculated from Q. In this case, the association contribution 
to the fugacity coefficient in (85) is non-zero so that the CPA EOS yields different results from 
the Peng-Robinson EOS on which it is based. The molar enthalpy of water as an ideal gas can 
be calculated using Q, with a = 3.470, b = 1.450 x 10“^ and c = 0.121 x 10® [Tl|. An 

alternative and slightly more accurate correlation for has been presented by Cooper US). Using 
the same reference point mentioned in the previous section (internal energy of liquid water along 
the VLE curve at 273.16 K), we calculate the enthalpy according to 


h{T, P) = -RT'^ 


/din 4>{T, P) 




dT 


+ h}^{T) - h{T = 273.16 K, P = 0.006117 bar), (4) 


P,n 


where 0.006117 bar is the vapor pressure of water at 273.16 K. Equation (|^ should also contain 
a term that represents the difference between the enthalpy and the internal energy of water at 
(T = 273.16 K, P = 0.006117 bar), but this term is so small that we neglect it. 

Results for five isotherms are shown in Eigures l^a) and [^b), while the molar enthalpy of 
vaporization is shown in Eigurej^c). Similar to the case for CO 2 , the the EOS parameters for water 


(such as e and k described in Section|B.2.3 1 are fit to vapor pressure data. Predictions of A/ivap are 


accurate to within a few percent as long as the temperature is not close to the critical temperature. 
Spycher and Pruess have obtained good agreement with experimental data for pure water, except for 
the enthalpy of the vapor along certain regions of the VLE curve [1|. A comparison of their model’s 
predictions with those of the CPA EOS along the VLE curve is shown in Eigure]^ The accuracy of 
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Figure 2: CPA EOS-predicted molar enthalpies of pure CO 2 (solid lines) vs. temperature along the 
vapor-liquid equilibrium (VLE) curve. Experimental data (dots) are obtained from NIST jl5| . The 
dashed vertical line crosses the abscissa at 304.14 K, which is the critical temperature of CO 2 . (a) 
Vapor-phase and liquid-phase enthalpies along the VLE curve, (b) presents the molar enthalpy of 
vaporization A/ivapj which is the difference between the vapor and liquid enthalpies in (a). 
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Eigure 3: CPA EOS-predicted enthalpies of pure H 2 O vs. pressure along five isotherms. Theoretical 
calculations are depicted in solid lines and are calculated from (|^. Experimental data (dots) are 
obtained from NIST m- The absolute error of the predictions in (a) is shown in (b). The 
1.0 kJ/mole error for the 500 K isotherm is associated with the vapor formed during the phase 
change at that temperature (see Figure Q. 
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the two models is similar in the range of temperatures between 373-550 K (less than 2.5 % maximum 
error for both). Above 550 K, the CPA EOS yields noticeably better results, although it fails to 
reproduce the sharp decrease in the vapor-phase enthalpy near the critical temperature. More 
importantly, however, Spycher and Pruess state that their model cannot calculate the enthalpy of 
water below 373 K. This could have practical implications for the applications described in Sectionj^ 
For example, in numerical simulations of CO 2 sequestration leakage, Spycher and Pruess’ model 
cannot be used to calculate the enthalpy in the aqueous environments of the shallow aquifers and 
nearby leakage pathways. 

3 Enthalpy of mixtures 

Mixtures of CO 2 and water can in general exist as two different phases, a C02-rich phase and an 
aqueous (H20-rich) phase. We treat the C02-rich phase as being composed only of CO 2 and water, 
while the aqueous phase may have dissolved sodium chloride. 

3.1 C02-rich phase (CO2 + H2O) 

At high (T, P) conditions, the C02-rich phase may be composed almost entirely of CO 2 that exists 
in a relatively dense supercritical state. The solubility of water in this form of CO 2 is less than 
2 mole % |2], and often times it is even much less than that. Since the mixture is almost pure 
CO 2 , the Peng-Robinson EOS is a very good approximation to the CPA EOS, and we have already 
demonstrated in Section [2.1| that it can accurately model the behavior of CO 2 over a wide range 
of (T, P), including supercritical conditions. However, at sufficiently low pressures for a given 
temperature, CO 2 may exist as a gaseous substance. When water vapor is exposed to CO 2 at 
these conditions, the resulting gaseous mixture may have very high concentrations of water, since 
gases are completely miscible in each other. These mixtures find applications in technologies such 
as oxy-fuel combustion, which is a method for carbon capture involving the formation of mixtures 
rich in CO 2 and water vapor [Hi- 

Experimental enthalpy data for CO 2 -H 2 O mixtures, which are much less common than data for 
the pure components, have been presented by Patel and Eubank m- They have studied C 02-rich 
gaseous mixtures between 323 K and 498 K where the water vapor mole fraction can be as high as 
50 %. The same set of data is also considered by Spycher and Pruess [T]. Patel and Eubank choose 
the reference to be the enthalpy of the corresponding ideal gas mixture at T = 273.16 K. With this 
reference, we calculate the molar enthalpy h{T, P, z) of the real fluid mixture from the CPA EOS 
according to 


h{T,P,z) = -RT^Ylzi ^ + /lig(T) - pg(r = 273.16 K), 

where h^T) = Zi is the mole fraction of component i, and the pure component 

ideal gas enthalpies h^f{T) are calculated as described in Section Enthalpy results for several 
isotherms of CO 2 -H 2 O mixtures at two different compositions are shown in Figure Spycher 
and Pruess’ model is more accurate at low pressures. As the pressure increases, the CPA EOS 
becomes more accurate. The point at which this switch occurs depends on the temperature and 
the composition; the CPA EOS tends to be more accurate at higher concentrations of water and 
at lower temperatures. In Figure |^a), the maximum error for the CPA EOS is 0.26 kJ/mole. In 
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Temperature (K) 

Figure 4: CPA EOS-predicted molar enthalpies of pure H 2 O vs. temperature along the VLE curve. 
Results from Spycher and Pruess [I] are included for comparison. Experimental data (dots) are 
obtained from NIST m- The dashed line vertical crosses the abscissa at 647.1 K, which is the 
critical temperature of H 2 O. A magnified view of the vapor-phase results in (a) is shown in (b). 
Predictions of the molar enthalpy of vaporization is compared with experimental data in (c). 
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Figure 5: CPA EOS-predicted molar enthalpies of gaseous CO 2 -H 2 O mixtures (solid lines) along 
different isotherms vs. pressure. Experimental data (dots) are obtained from Patel and Eubank |18] . 
Results from Spycher and Pruess (dashed lines) [I] are included for comparison, (a) 9:1 mole ratio 
of CO 2 to H 2 O; (b) 1:1 mole ratio of CO 2 to H 2 O. 


Eigure|^b), the maximum error is 0.17 kJ/mole. Spycher and Pruess state that their maximum 
error is about 0.40 kJ/mole (13 kJ/kg). Besides the fact that the CPA EOS can simultaneously 
predict phase equilibria behavior, density, and enthalpy, its main advantage is that it can calculate 
the enthalpy of pure water or mixtures involving water below 373 K. Although not shown in 
Eigurej^a), Patel and Eubank present a few enthalpy data points at 348 K for the 90 % CO 2 
mixture. The maximum error of the CPA EOS for these points is less than 0.03 kJ/mole. 


3.2 Aqueous phase (H2O + CO2 + NaCl) 

To the best of our knowledge, there are no reported molar enthalpy measurements for aqueous 
mixtures with CO 2 . Measurements of the CO 2 molar enthalpy of solution A/igoi are available, 
however. This quantity is defined as 


A/isoi = Hco2{T,P,zco2 -t 0) - hco2{T,P). 


(5) 


That is, A/isoi is the difference between the partial molar enthalpy of CO 2 in an infinitely dilute 
aqueous solvent and the molar enthalpy of pure CO 2 . Koschei et al. [H] have presented Ahgoi 
data at 323 K and 373 K in three different aqueous solvents: pure water, 1 molal (m) NaCl, 
and 3 m NaCl. As explained in Section B.2[ the CPA EOS is limited to aqueous-phase mixtures 
containing only pure water. Eor such mixtures, Eigurej^ compares Ah^oi calculations from the CPA 
EOS with experimental data from Koschei et al. The agreement at 323 K is poor, while that at 
373 K is relatively good. Since the EOS accurately calculates the molar enthalpy of pure CO 2 , 
most of the error in Ahso\ comes from calculation of the infinite dilution partial molar enthalpy. 

Despite the large error in A/igoi, the CPA EOS can still provide accurate predictions of the 


molar enthalpy h{T,P,z) of the mixture. This can be reasoned as follows. From (26), the molar 
enthalpy of H 2 O-CO 2 mixtures is 










Figure 6: (a) CPA EOS-predicted CO 2 molar enthalpy of solution in pure water vs. pressure at 
two different temperatures. Experimental data are obtained from Koschei et al. |19] . (b) shows 
the absolute error of the predictions. Despite the large error in A/igoi at 323 K, the EOS can still 
accurately compute the molar enthalpy h{T, P,z) of H 2 O-CO 2 mixtures, as explained in the text. 


h{T,P,z) = ZCO 2 HCO 2 + ^HaO^^HzO- 

Eor the temperatures and pressures considered in this report (say, T < 500 K and P < 500 bar), 
the CO 2 solubility in pure water is less than 3.5 mole % [2]. Suppose that the error in the cal¬ 
culation of the partial molar enthalpy Hqq^ is 4.0 kJ/mole CO 2 , which is a value greater than 
the maximum error of 3.5 kJ/mole CO 2 in Eigure Then the error in the product zco 2 ^C 02 ^ 
assuming that ^002 = 0.035, is 0.14 kJ/mole. The partial molar enthalpy of water is similar 

to the molar enthalpy h-n^o of the pure component because the aqueous mixture is predominantly 
water. Eor example at 323 K, the CPA EOS predicts that the difference between the two quantities 
is less than 0.006 kJ/mole H 2 O even for C02-saturated water at 500 bar. It is therefore reasonable 
to expect the error in .ffH 2 0 to be similar to that for /iH 2 0 shown in Eigurej^ where the maximum 
error is 0.25 kJ/mole H 2 O in the liquid phase. Thus, the error in h(T,P,z) over a wide range of 
conditions will likely be less than 0.40 kJ/mole (~ 0.14 kJ/mole -|- 0.25 kJ/mole), which is rather 
small considering the expected magnitude of h(T,P,z). Nonetheless, we propose a generalized 
version of the Li and Eiroozabadi CPA EOS in Section and explain how it may provide better 
agreement with the excess enthalpy, which is related to Akg^i. Another alternative is the Duan and 
Sun CO 2 activity coefficient model [H], which not only agrees better with A/igoi measurements in 
pure water, but also provides an accurate way to calculate the CO 2 partial molar enthalpy (and 
subsequently A/igoi) in brine solutions. Their study is the focus of the rest of this section. 

Duan and Sun developed their model for the purpose of performing CO 2 solubility calculations 
in pure water and brines. Appendix 0 describes how it can also be used to calculate the CO 2 
partial molar enthalpy Hco^{T,P,vn) and A/igoi. Here, m represents the molalities of all the 
solutes dissolved in the aqueous solvents. Eigure compares A/isoi predictions from their model 
with data from Koschei et al. Eor mixtures in pure water (m = 0), there is marked improvement 
over the CPA EOS at 323 K. Since CO 2 is even less soluble in brines than in pure water, we follow 
the same reasoning in the preceding paragraph to argue that the CO 2 contribution to the error 
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Figure 7: CO 2 molar enthalpy of solution derived from Duan and Sun’s CO 2 activity coefficient 
model [3] vs. pressure at 323 K (a,b) and 373 K (c,d) and three different NaCl molalities m. 
Experimental data are obtained from Koschei et al. [Kl¬ 


in the molar enthalpy h(T, P, m) of an H20-C02-NaCl aqueous mixture is relatively small as long 
as Hco^iT^P^-m) is predicted to within a few kJ/mole CO 2 . The accuracy of Duan and Sun’s 
model is well within this limit, at least for the conditions where experimental data are available 
(temperatures, pressures, and NaCl concentrations up to 373 K, 200 bar, and 3 molal, respectively). 

In order to obtain the molar enthalpy h(T, P, m), we must calculate the contributions from water 
and NaCl. Rather than treating water and NaCl separately, one approach is to treat them together 
as a single pseudocomponent, which we label as brine. In this approach, the H20-C02-NaCl mixture 
may be thought of as a binary brine-C02 mixture. A correlation for the specific enthalpy /ibrine of 
sodium chloride solutions, in units of kJ/kg brine, has been presented by Michaelides [2D]. This 
correlation contains an error in one of the coefficients that was corrected by a later study |2I] . For 
every kilogram of water in the brine-C02 solution, the total moles are 

1000 

— -h mNaCl + "^002 5 

J0'H20 
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and the total mass of brine (H 2 O + NaCl), in units of kg brine/kg H 2 O, is 


1 + JT^NaCl 


-/i^NaCl 
1000 ^ 


where MNaCi = 58.44 g/mole is the molecular weight of the salt. As a result, the enthalpy of the 
solution per kilogram of water is approximately given by 


mco 2 Hco 2 {T, P, m) + + mNaCl 

The molar enthalpy is therefore 


AfNaCl ^ 

1000 J 


hhv 


h(r,P,m) 


^002 

1000 

vy-1- "TNaCl + mc02 

L MH 2 O 


Hc02iT,P,ui) + 


^1 + mNaCl 


AfNaCl \ 
1000 ) 


TOOO' 

yy-1- m-NaCl + "^002 

L IWH 2 O 


^brine- (6) 


Equation Q is an approximation because h{T, P, m) should involve the specific enthalpy of brine 
that contains dissolved CO 2 , and not the specific enthalpy of brine by itself. However, since the 
mixture is predominantly brine, the two quantities are not expected to be very different, just like 
how HhjO and / 1 H 2 O are similar. Michaelides states /ibrine is accurate to within 3 %. This value 
is comparable to the error of the liquid water molar enthalpy in Figure Since Duan and Sun’s 
model can compute Hqq^ more accurately than the CPA EOS for most of the conditions, it is 
reasonable to expect h{T, P, m) predictions to deviate from experimental data by no more than 
0.40 kJ/kg (the maximum estimated error for the CPA EOS) over a wide range of conditions. 


4 Excess enthalpy 


This section explains how relaxing two of the assumptions in Li and Firoozabadi’s version of the 
CPA EOS may allow for better agreement with experimental data for the molar excess enthalpy. 
This quantity is distinct from the molar enthalpy that we have so far considered. From the general 


definition in (24), the molar excess enthalpy is 


^excess^y^ P, z) = h{T, P, z) - P^{T, P, z), 

where h^^(T, P, z) is the molar enthalpy of the corresponding ideal mixture at the conditions 
(T,P,z). Since the molar enthalpy of an ideal mixture is calculated from the pure component 
enthalpies (see Section A.3|, A/isoi in (§ is equivalent to the partial molar excess enthalpy of CO 2 
in an infinitely dilute mixture. Therefore, the generalized CPA EOS proposed in this section is 
expected to also improve the agreement depicted in Figure Bottini and Saville have measured 
the molar excess enthalpy of gaseous CO 2 -H 2 O mixtures in the temperature range 520 - 620 K, 
pressure range 14 - 45 bar, and CO 2 mole fractions between 20 - 78 % [52]. For these mixtures, 
the CPA EOS yields an average absolute error of 1.1 kJ/kg (this corresponds to a 13.3 % error) 
and a maximum error of 5.6 kJ/kg. Figure compares CPA predictions of the molar excess en¬ 
thalpy of gaseous CO 2 -H 2 O mixtures with data from two other studies. The absolute error can be 
a significant fraction of the excess enthalpy, especially at the higher temperatures. Nevertheless, 
these results are an improvement over those from Spycher and Pruess, who state that their average 
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Figure 8: CPA EOS-predicted molar excess enthalpies of gaseous CO 2 -H 2 O mixtures vs. pressure 
at three different T. The experimental data at 598 K are from |23j and represent 1:1 mole ratio 
mixtures of CO 2 to H 2 O. Data at 694 K and 803 K are from [25], where the composition ranges 
between 47-51 mole % of C02. (b) shows the absolute error of the results in (a). 


absolute error for 1:1 gaseous mixtures is about 0.30 kJ/mole (10 kJ/kg), and their average error for 
the Bottini and Saville study is 2.4 kJ/kg [I]. We note that an earlier study of H 2 O-CO 2 mixtures 
has achieved good agreement with excess enthalpy data up to 598 K, although their model also 
exhibits significant errors at higher temperatures |23|. This study uses optimization techniques to 
fit parameters to many different sets of data. The generalized corresponding states principle is 
employed to make the parameter space more manageable. Their model is targeted towards only 
the aqueous phase, however. It is not applicable to pure CO 2 or to C02-rich mixtures with water 
vapor, like the 90 mole % CO 2 mixtures in Figure |^a). 

We may explain the behavior in Figure by considering the nature of intermolecular forces, 
particularly hydrogen bonding, that exist in mixtures of CO 2 and H 2 O (see Sections B.l and B.2|. 
Breaking hydrogen bonds requires an energy input on the order of 10 kJ/mole [25]) which is at least 
an order of magnitude larger than the excess enthalpy values. The large disparity suggests that the 
excess enthalpies are rather small in magnitude, so that even small discrepancies between theory and 
experiment appear as large relative errors in Figure |^b). The discrepancy may be largely due to the 
two assumptions listed in Section B.5, The assumptions state that the a (hydrogen-bond donor) and 
(3 (hydrogen-bond acceptor) sites are equally numerous on a molecule of component i and associate 
symmetrically so that their site fractions Xia and Xip can be described by a single function Xi- In 
order to see the consequences of these assumptions, let us examine, as a representative example, a 
1:1 CO 2 -H 2 O gaseous mixture at 598 K and 66.5 bar. The relevant values are: 


• Experimental molar excess enthalpy = 0.73 kJ/mole. 

• CPA-predicted molar excess enthalpy = h(T,P,z) — /i™(T, P, z) = 0.92 kJ/mole. 

• CPA-predicted molar enthalpy departure function = h(T, P, z)—P^(T, P, z) = —1.46 kJ/mole. 
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Qualitatively, the signs and relative magnitudes of the excess enthalpy and departure function are 
as expected. The departure function should be negative because the intermolecular interactions 
present in the real (non-ideal) mixture reduce the energy of the mixture compared to the non¬ 
interacting ideal gas mixture. The reason why the excess enthalpy is positive can be understood 
by considering the process of mixing CO 2 with water. Mixing involves three processes: 

1. Breaking CO 2 -CO 2 intermolecular forces, primarily London dispersion forces. 

2. Breaking H 2 O-H 2 O intermolecular forces, the most prominent being hydrogen bonds. 

3. Forming CO 2 -H 2 O intermolecular forces, the most prominent being hydrogen bonds. 


The first two require energy input while the last one releases energy. The CO 2 -H 2 O hydrogen bonds 
are not taken into account in the ideal mixture, since the ideal mixture enthalpy is calculated from 
the pure component enthalpies. These hydrogen bonds are weaker than H 2 O-H 2 O hydrogen bonds 
because C is more electronegative than H, so the lone pairs of electrons on the oxygens in CO 2 are 
less effective hydrogen-bond acceptors than the lone pairs on the oxygen in H 2 O. As a result, the 
real mixture should be higher in energy than the ideal mixture, meaning a positive excess enthalpy. 

A molecule of CO 2 has four /3 sites, representing the 4 lone pairs on the two oxygens, and no 
a sites (hydrogens). Treating the a and /3 sites as symmetric and equally numerous per molecule 
will overpredict the number of CO 2 -H 2 O hydrogen bonds formed since CO 2 does not have any 
a sites. Consequently, overprediction of the number of CO 2 -H 2 O hydrogen bonds will lead to an 
underprediction of the number of H 2 O-H 2 O hydrogen bonds in the mixture. This is because if a 
H 2 O molecule forms a hydrogen bond with CO 2 , it has fewer sites available to form a hydrogen bond 
with another H 2 O. As a result, the theoretically-predicted excess enthalpy will be more positive 
than the experimental value, as we have observed. The small magnitude of the excess enthalpy 
is also consistent with our discussion above. We have reasoned that is positive because 

CO 2 -H 2 O hydrogen bonds are weaker than H 2 O-H 2 O hydrogen bonds. This is due to the fact 
that carbon is more electronegative than hydrogen. However, since carbon is only slightly more 
electronegative than hydrogen, we would expect the excess enthalpy to be quite small, certainly 
much less than the energy required to break a hydrogen bond, which is 0(10 kJ/mole). 

In order to achieve better agreement with excess enthalpy measurements, Li and Firoozabadi’s 
CPA EOS must be modified so that it distinguishes between a and /3 sites, treating them asymmet¬ 
rically. We suggest relaxing the two assumptions in Section |B.5| so that equations in Sections |B.3 


and B.4 for the fugacity coefficients (l)i (on which all EOS predictions of the density, solubility, and 
enthalpy are based), pressure, and compressibility factor Z are solved instead. These equations 
are derived from the departure function in (60) or equivalently, the one in (51). From a practical 
perspective, the asymmetric treatment will introduce additional cross-association parameters (see 
Section B.2.3) that can be fit to excess enthalpy data. The pure component parameters that exist 
in the current model, such as e and k described in Section IB. 2. 31 do not need to be modified. 


5 Conclusions and future work 

We have shown that the CPA EOS accurately calculates the molar enthalpy of pure CO 2 , pure 
water, C02-rich mixtures, and aqueous mixtures for the conditions relevant to the applications 
discussed in Section Since the CPA EOS has previously been demonstrated to provide good 
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agreement with density and solubility data [2], it can therefore serve as a single, unified model for 
thermodynamic property predictions of all of these fluids. Compared to the Spycher and Pruess 
model [T], the CPA EOS more accurately calculates the vapor-phase enthalpies of pure water at 
temperatures above 550 K along the vapor-liquid equilibrium curve (Figure [^. For CO 2 -H 2 O 
mixtures, the CPA is more accurate at lower temperatures, higher pressures, and higher water 
concentrations (although both models predict the molar enthalpy data in Figure well). For 
practical purposes, the most important difference between the two models is that the CPA EOS is 
not limited only to temperatures above 373 K for pure water or water-containing mixtures. 

The CPA EOS is not applicable to aqueous H20-C02-NaCl mixtures, however. For these 
brine mixtures, we suggest calculating the molar enthalpy with We have argued that is 
expected to be accurate to within 0.40 kJ/mole, at least for the conditions where enthalpy data 
are available jl9| (temperatures, pressures, and NaCl concentrations up to 373 K, 200 bar, and 
3 molal, respectively). Equation ^ involves the partial molar enthalpy of CO 2 , which we obtain 
by evaluating temperature derivatives of Duan and Sun’s CO 2 activity coefficient model [3]. In a 
similar fashion, one can obtain the partial molar volume of CO 2 by evaluating pressure derivatives. 
The density of the brine can then be calculated from this partial molar volume. However, we find 
that the CO 2 partial molar volume derived in this way does not closely agree with experimental 
data m- Density calculations require finer precision than enthalpy calculations because the density 
increase from CO 2 dissolution into water or brine is small. It is less than 2 % even at pressures where 
the CO 2 solubility is high (3 mole %). As a result, we recommend the density of H20-C02-NaCl 
mixtures be calculated using the models presented in more recent studies |28l 123] . 

The CPA EOS may be combined with Duan and Sun’s model to form a 'y-4> model for H 2 O- 
C02-NaCl mixtures. In this approach, the CPA EOS is used to compute the density, enthalpy, 
and composition of the C 02 -rich phase, which we treat as being composed only of CO 2 and water. 
For the aqueous phase, Duan and Sun’s model is used to calculate the CO 2 solubility and the CO 2 
partial molar enthalpy. The partial molar enthalpy may be combined with the specific enthalpy of 
brine (computed from the studies |2UII21] described in Section 3.2) to obtain the molar enthalpy (§ 
of the mixture as whole. The density of the aqueous phase may be calculated from one of the 
aforementioned models [2E1I23]. A more theoretically robust alternative to a ’’y-cj) model is a (ji-cj) 
model, in which the same EOS is used for both phases. In our future work, we plan to develop a 
model by modifying the CPA EOS so that it can accurately predict the thermodynamic properties 
of both phases. The first step is to generalize the CPA EOS as described in Section]^ to further 
improve the agreement with CO 2 -H 2 O mixture enthalpy data. After this step, we may extend the 
applicability of the CPA EOS to H20-C02-NaCl mixtures by adding more terms to the departure 
function (32) to account for intermolecular forces such as ion-ion and ion-dipole interactions. This 
work will involve many challenges, some of which arise due to the lack of experimental molar 
enthalpy data for aqueous mixtures. Nonetheless, a systematic and rigorous approach that builds 
on previous studies of electrolytic solutions [Ml ED [32] may prove to be successful. 
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A Basic concepts in the thermodynamics of fluid mixtures 

A.l Chemical potential, fugacity, fugacity coefficient 

Many fundamental thermodynamic relations for real fluid mixtures are inspired by the correspond¬ 
ing relations for ideal gas mixtures. One such example is the relation between the chemical potential 
P, n) of component i in an ideal gas mixture and the pure component ideal gas chemical po¬ 
tential Here, T is the temperature, P is the pressure, n represents the mole numbers of 

the mixture’s components. The difference between P, n) and P) is 


/i‘®(r, P, n) - /if (T, P) = RT In Zi, (7) 

where R is the gas constant, and Zi = rii/N is the mole fraction of component i. Equation 0 can 
be derived formally using statistical thermodynamics, but it can also be derived in a more informal 
manner as follows [33] . Suppose a pure component exists as an ideal gas at pressure P. If this pure 
component is mixed isobarically with other ideal gases to form a mixture whose pressure is also P, 
the pressure of the component in the mixture becomes the partial pressure ZiP so that 


/if(r,p,n)-/if(r,p) 




T,n 


dP'. 


In our study, we consider only systems where external fields like gravity are negligible and there is 
only pressure-volume work. For such systems. 


C 

dG = -SdT + VdP + '^Hidni, ( 8 ) 

i=l 

where G is the Gibbs energy, S is the entropy, V is the volume, and c is the number of components. 
Applying the equality of mixed partial derivatives to ([^, we have 


rdv 

\dni 


T,P,n, 



(9) 


where the subscript rij in the derivative of V means that all mole numbers besides rii are held fixed. 
Using ([^ and the volume U’® = NRTjP of an ideal gas, we obtain 






dui 


T,P',n, 


dP' = RT 



RT in Zi , 


which leads to Q. Similarly, the chemical potential /if(T, P, n) is related to /if(T, P',n) at a 
different pressure P' (but same temperature and composition) by 




dP” = 


T,n 


fif) 

Ip' y orii J 


dP” = RT In 


T,P",n, 


P 

Y' 


( 10 ) 


Combining ([(^ and (101, the difference between /rf (T, P, n) and P', n'), which is at a different 

pressure and composition but same temperature, is 


/xf (T, P, n) - /xf (T, P', n') = [/xf (T, P) - /xf (T, P')] + RTln = RTln ^. (11) 
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In analogy to (0, for any fluid (not just ideal gases) each component i is assigned a fugacity 
fi = fi{T,P,n) so that the difference between ^i{T, P, n) and P', n') is 


fii{T,P,n) - fii{T,P',n) = RTIn 


n{T,P,n) 

MT,P',n')- 


( 12 ) 


In order for (12| to be consistent with (111, the fugacity of component i in an ideal gas mixture 
must be equal to its partial pressure 


fl%T,P,n) = ZiP (13) 

The difference in the chemical potential of component i in a real fluid and the chemical potential 
of that same component in an ideal gas mixture at the same (T, P, n) is 

/ii(r,P,n) -/rf(r,P,n) = RTln = RTlncj^i, (14) 

where 

0,(r,P,n)= (15) 

ZiP 

is the fugacity coefficient of component i. Since any fluid mixture behaves like its corresponding 
ideal gas mixture as the pressure P decreases to zero, fi becomes P, n) = ZiP and (/>* becomes 

unity in this limit: 


lim fi(T, P, n) 


ZiP, 


hm (()i(r,P,n) = 1. 
P—>0 


If the fluid in question is a pure component fluid. 


lim f(T, P) = P, 

lim d>(T, P) = 1. 

P-s-O 

For a multiphase system of p phases to be in equilibrium, a fundamental thermodynamic require¬ 
ment j34j is the equality of chemical potentials 


12 P 

= Mi =■■■ = Pi 


for all i. 


(16) 


Using (12), the condition for phase equilibrium in (16) can be equivalently stated as the equality 
of fugacities 


fi = fi =■■■ = fi for all i. 


(17) 


The fugacity fi can be easily calculated from the fugacity coefficient (^i and vice-versa using (15). 
If we know how the fugacity of each component varies as a function of temperature, pressure, and 
composition, (0 can be solved to determine the equilibrium composition (i.e., the solubility of each 
component) in each phase. Furthermore, we will see in Section (A.31 that other thermodynamic 
functions, like the density and the enthalpy, can in principle be derived from the fugacities or the 
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fugacity coefficients. In essence, the purpose of an equation of state is to provide a relation from 
which we can derive the functional form of the fugacities for all components. Consequently, the 
equation of state can be used to derive other thermodynamic functions of the mixture as well. In 
this sense, the equation of state provides a complete thermodynamic specification of the fluids to 
which it is applicable. 


A .2 Activity coefficient, activity, excess property, partial molar property 

In analogy to Q, an ideal mixture (also called an ideal solution) is defined to be one in which 

= RTlnzi. (18) 


Here, iJ,i{T,P) is the chemical potential of pure component i in the real fluid state, as opposed to 
the hypothetical ideal gas state in ([^. From (12) and (181, the fugacity /)™(r, P, n) in an ideal 
mixture and the fugacity fi{T, P) of the pure component are related by 


friT,P,n) = zJiiT,P). (19) 

An ideal gas mixture can be thought of as an ideal mixture where the pure component fugacity 
fi{T, P) = P for all conditions such that the fugacity fl^{T, P, n) in the mixture equals the partial 
pressure ZiP. In other words, an ideal gas mixture is an ideal mixture where (13 1 is satisfied. 
Equation (19) is a simple relation between a mixture property /™(r, P, n) and a pure component 
property fi(T,P)-, it is a key feature of ideal mixtures. A similar relation can be formed for a real 
(non-ideal) fluid mixture by introducing an activity coefficient 'ji = 7 j(T, P, n) defined such that 


MT,P,n) = zai{T,P,n)fi{T,P), 


( 20 ) 


where 


lim 7 j = 1 . 

Zi —>-1 

The activity a* of component i is defined as 


_ /.(r,p,n) 
a. 7^z^ , 

so that for any fluid. 


( 21 ) 


( 22 ) 


//i(r,P, n) -/ij(r,P) = PT In a*. (23) 

With these definitions, one can view an ideal mixture as a real fluid in which 7 * = 1 (or equiva¬ 
lently, aj = Zi) for all components over all conditions. An excess property is defined as 


pexcess^j.^ P, n) = E{T, P, n) - P™(r, P, n), (24) 

where E is an extensive property (e.g., Gibbs energy) of the mixture. 

A related concept to the excess property is the partial molar property. The partial molar 
property of E with respect to component i is denoted as Pj, and is defined as 


Ei = 



(25) 
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One example that we have seen previously in 


is the partial molar volume 


= 


_ /dV 

[dui 


T,P,n, 


Another example is the chemical potential in (l8|), which is the partial molar Gibbs energy 


— Gi 



There exists a very important relation between E and its partial molar properties Ei. It can be 
shown [321 ESI ES] that 


E{T,P,n) = '^mEi, 

i=l 

or on a molar basis, if e(r, P, z) = E{T^ P, n) /N^ 

c 

e{T,P,z) ='^ZiEi. (26) 

i=l 

We denote the set of all mole fractions as z = n/N. The molar Gibbs energy is 


9{T, P,z) = ^iGi = XI 

i=l 


(27) 


i=l 


Due to complexities arising from intermolecular interactions, the partial molar property of a mixture 
is in general not equal to the molar property ei{T,P) = Ei(T, P)/ni of pure component i. That 
is, Ei / ei{T,P) except i n spe cial cases. One such special case is the partial molar enthalpy of 

^ \T,P,n) = hi{T,P). We will 


an ideal mixture. Section 


A.3 


shows that for an ideal mixture, i7, 


im/ 

i 


later use this relation to calculate the enthalpy of real fluid mixtures. It is important to note that 
(T, P, rij) must be held fixed in the derivative for it to be considered a partial molar property. For 

that fii = {dP/dni)\j,y^^,, 


B.3 


where F is the Helmholtz energy of 


example, we will see in Section 
a fluid. Since it is the volume, and not the pressure, that is held fixed we conclude that Hi ^ Ei 
Instead, the partial molar Helmholtz energy is defined as 


Fi = 


dF 

drii 


T,P,ni 


We can also define partial molar properties of excess properties. For example, 

/ Q^exce, 


Gf 


\ drii 


T,P,ni 


= /r*(r,P,n)-/xnr,H,n). 


Using this relation along with (12) and (19|“(20), we have 


^excess ^ p^ ^ 


fi{T,P,n) 

fr{T,P,n) 


RT In 


ZilifiiT, P) 

' Zifi{T,P) ' 


RT In 7 j. 


This result shows that if we have a smooth function (such as a polynomial) that describes the 
variation of the excess Gibbs energy with the composition, we can take composition derivatives 
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(i.e., evaluate the partial molar excess Gibbs energy to obtain the activity coefficient yj. In 

turn, we can use 7 * to compute the real fluid mixture fugacity fi{T, P, n) from the pure component 


fugacity fi(T,P) using (20). The smooth function is called the activity coefficient model. It is 
essentially a correlation which provides a continuous fit to a discrete set of excess Gibbs energy 
data, which are obtained through experimental measurements at different compositions. See j34ll37| 
for extensive discussions on activity coefficient models for various types of mixtures. Duan and Sun 
have presented a model for the GO 2 activity coefficient in mixtures where GO 2 is dissolved in 
aqueous sodium chloride solutions [3] . This model is partly based on the activity coefficient model 
of Pitzer, which is widely used to compute the excess Gibbs energy of electrolytic solutions 
The Duan and Sun model is the focus of Section 3^ and Appendix 


A.3 Enthalpy and density 

We stated in Section A.l| that the enthalpy of a fluid mixture can be derived from the fugacity 
coefficients (pi. To show this, we use the relation G = H — TS between the Gibbs energy G and 
the enthalpy H so that the partial molar properties are related by 


From ®, we have 


which leads to 


This result can be rearranged to 


Gi = ^ii = Hi- TS^. 


Si = 


dS 

dm 


T,p,n, 


/J-i — FIj + T 


dm 

dT 


d / m 
dT \T 


dm 

dT 


P,n 


P,n 


Hi 

J'2 ' 


P,n 


(28) 


Since an ideal mixture is defined by (181 



_ d 

f m{T,P) + RTIn Zi\ 

_ d 


P,n 

1 T ) 

P,n dT 

1 T ) 


P,n 


This is the result we alluded to in Section 


A.2 


the partial molar enthalpy FF, 


hiiT,P) 

T2 ■ 

(29) 

is equal to the 


molar enthalpy hi{T,P) of the pure component. The enthalpy of ideal mixtures, which includes 
ideal gas mixtures, can thus be computed from the pure component enthalpies. In other words, 
ideal mixtures are characterized by zero enthalpies of mixing: 


AhZ. = F-(T, P, z) - ^ zMT, P) = E - h, 


i=l 


2=1 


= 0 . 


Applying (28) and (29) to (14), we have 
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A ( Mi - 
dT \ T 


P,n 


^ fdlnUT,P,n) 

= - df - 


H,-Hf _ H,-h^ 


P,n 


'J'2 


rp2 


Combining this result with (26), the molar enthalpy of any fluid mixture is 

( 91n(/)j(r,P, n) 


h{T, P,z) = J2 ZiHi = -RT^ E 1^- 


dT 


+ E 

i=l 


(30) 


2=1 2=1 

This is the relation we seek between the fugacity coefficients and the enthalpy of the mixture. We 
can calculate the enthalpies /i)® of the pure components in their ideal gas states using correlations 
found in the references specified in Section 

The density can be derived from the fugacity coefficients as well. Using along with ( |14[ ) 
and (26), the molar volume v = V/N is computed from the partial molar volumes U according to 


C C 

U = E = E 


f dm 
* [dP 


2 = 1 2 = 1 

For an ideal gas where = NRT/P 


Tn j=i 


dP 


I T>rr f din (j)i 


dP 


= = 
' -j 




r,n 


dui I 




RT 

IP' 


Therefore, 


and the molar density p is 


= RTY,Zi 


2=1 


1 fd\n(t)i 
\ dP 


T,n 


P = - = 
V 


RTE^=lZ^ 


1 fd\n(t)i\ 
P^ \ dP ) 


T,n 


The mass density is 


Pmass — P 'y 


2=1 


where Mi is the molecular weight of component i. 


r,n 


B Cubic-plus-association (CPA) equation of state 

B.l Intermolecular forces and departure functions 


It was stated at the end of Section A.l that the fundamental role of an equation of state (EOS) is to 
provide an expression from which all other thermodynamic functions of a fluid can be derived. For 
equations of state that model complicated fluids like mixtures of CO 2 and water, this expression is 
often in the form of a departure function [33 E3], which is sometimes also referred to as a residual 
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or a residual function (HU ES] • A departure function is defined as the difference in some property 
(e.g., Helmholtz energy) of a real fluid and that same property of the fluid if it were to exist as an 
ideal gas mixture. This definition resembles the definition of an excess property in (24), but there 
are two important differences. First, the reference fluid in (24| is an ideal mixture, while it is an 
ideal gas mixture in a departure function. Second, excess properties are always defined as differences 
in some property E for two different fluids at the same temperature, pressure, and composition 
(T,P,n). In contrast, departure functions can involve two different fluids at the same (T, P, n), 
or two fluids at the same (T, y,n). The two types of departure functions are not necessarily the 
same. An ideal gas is by definition composed of non-interacting molecules represented by volumeless 
entities (points). In a real fluid, the departure from ideality arises because the molecules have a 
non-zero volume and interact with each other through intermolecular forces. The volume of the 
molecules is itself a reflection of repulsive intermolecular forces that have their origins in electrostatic 
forces and quantum mechanical effects (i.e., the Pauli exclusion principle) j39j. 

In order to understand the physical behavior represented by the CPA EOS, and indeed equations 
of state in general, we briefly review some basic concepts regarding intermolecular forces. The 
discussion here will also motivate the future work described in Sections |4] and [H Intermolecular 
interactions can generally be divided into two sets: physical and chemical/quasi-chemical. Physical 
interactions include London dispersion forces (which can be thought of as instantaneous dipole- 
induced dipole interactions), dipole-dipole interactions (also called Keesom forces), and dipole- 
induced dipole interactions (Debye forces). They are sometimes collectively referred to as van der 
Waals forces [SB]. On a per-mole basis, the van der Waals forces are weaker than other types of 
interactions. However, because they form rather easily, especially London dispersion forces, they 
are relatively numerous in a given fluid and can significantly influence the properties of a fluid. 

Chemical interactions include hydrogen bonds, which act between a molecule that contains the 
electronegative atoms oxygen (O), nitrogen (N), or fluorine (F) bound to hydrogen (H) and another 
molecule that contains O, N, or F. The hydrogen develops a partial positive charge, while the 
electronegative atom to which it is covalently bonded develops partial negative charge. Hydrogen 
bonding allows molecules of a fluid to associate to form polymers. For example, methanol consists of 
CH 3 OH monomers, some of which polymerize to form CH 3 OH dimers, trimers, or longer-chained 
polymers. The monomer units are linked together (associated) by hydrogen bonds between the 
oxygen on a methanol molecule and the hydrogen in the hydroxyl group of another methanol. 
Hydrogen bonds are classified as chemical interactions because they are quite strong and resemble 
chemical (i.e., covalent) bonds in that they involve partial overlap of the electron clouds of the 
atoms involved in the interaction [M]. In contrast, physical interactions do not involve significant 
overlap of electron clouds. 

Hydrogen bonding plays an essential role in determining water’s unusual properties. Ice consists 
of a regular lattice of water molecules, with each molecule being bound to four other molecules 
through hydrogen bonds. Each molecule is said to have four associating sites (Figure]^, which 
come in two varieties: 2 hydrogen-bond donor (a) sites representing the two hydrogens, and 2 
hydrogen-bond acceptor (/I) sites representing the two lone pairs of electrons on the oxygen. In 
liquid water, the molecules also form hydrogen bonds with each other, but these interactions are 
more transient; hydrogen bonds are constantly forming and breaking, and the molecules are able to 
move past each other. The regular, hydrogen-bonded lattice structure of ice limits the maximum 
packing achievable, which is why ice is less dense than liquid water at the normal melting point. 

In fluid mixtures, the components can be either self-associating or non-self-associating, based on 
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Figure 9: Hydrogen bonding in water. Each molecule can form up to four hydrogen bonds since 
it has four associating sites: 2 hydrogen-bond donor (also called electron acceptor) sites represent¬ 
ing the two hydrogens, and 2 hydrogen-bond acceptor (electron donor) sites representing the two 
electron lone pairs on the oxygen. We will refer to these two site types as a and /3 sites, respectively. 


whether the pure component forms hydrogen bonds. In mixtures with at least one self-associating 
component, there may be cross-association between the self-associating component and the other 
components. For example, CO 2 may cross-associate with water because the lone pairs on the 
oxygens of CO 2 can serve as hydrogen-bond acceptor sites. Molecules without associating sites, 
such as H 2 S or alkanes, may also cross-associate with water. In this case, the cross-association 
provides a simple way to implicitly account for the complicated physical and chemical processes 
that occur during the solvation of these molecules in water. 

B.2 Helmholtz energy departure function 
B.2.1 Overview of Peng-Robinson and CPA EOS 

The Peng-Robinson (PR) EOS belongs to the family of cubic equations of state, which are all based 
on the van der Waals equation of state [33] . The PR EOS is widely used in the oil and gas industry 
to model fluid mixtures containing only non-self-associating components, such as hydrocarbons. 
It, along with other cubic EOS variants like the Soave-Redlich-Kwong (SRK) EOS, involves an 
interaction parameter a, which accounts for van der Waals interactions, and a co-volume parameter 
b, which accounts for the non-zero volume of the fluid molecules. If a = 5 = 0, the PR EOS reduces 
to the ideal gas EOS. Since a accounts for only van der Waals interactions, the PR EOS in its 
original form cannot model fluid mixtures with associating components, such as water. 

One approach to model mixtures with associating components is to append the PR EOS, or 
a different cubic EOS, with an additional term that accounts for the association. This approach 
forms the basis of the cubic-plus-association (CPA) EOS. The CPA EOS in its modern form was 
introduced by Kontogeorgis et al. [ID] , although earlier authors have developed a similar EOS m- 
Other models for associating fluids besides CPA are described in various studies |l2l |l3| , and will 
not be considered further in this report. The cubic part of the CPA EOS by Kontogeorgis et al. 
comes from the SRK EOS. The association (hydrogen bonding) part is based on the thermodynamic 
perturbation theory of Wertheim mmsi SSI szj- In fluids without associating components, the 
CPA EOS of Kontogeorgis et al. reduces to the SRK EOS. 
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Li and Firoozabadi have presented a version of the CPA EOS where they introduce a novel 
scheme to model the cross-association between molecules of different components [2]. Their study 
concludes that cross-association is necessary to faithfully reproduce experimental vapor-liquid equi¬ 
libria (VLE) data, especially along the vapor phase, for mixtures containing water and one or more 
of the following components: CO 2 , H 2 S, and hydrocarbons (both aliphatic chains and aromat¬ 
ics). This conclusion is in agreement with earlier studies, who also report that cross-association 
is important in water/alkane mixtures and water/C 02 mixtures |481 112 ] , Li and Eiroozabadi’s 
original version of the CPA is targeted toward mixtures where water is the only self-associating 
component [2]. They later modified their EOS to model mixtures where asphaltenes serve as the 
self-associating component [49 [ I5U] . The Li and Eiroozabadi CPA EOS with water as the only 
self-associating component will be the focus of the present study. 

The starting point is to consider the CPA EOS’s expression for the departure function 
of the Helmholtz energy. Since the Helmholtz energy E is a natural function of (T, V, n), is 

defined as the difference between the Helmholtz energy of the real fluid and the Helmholtz energy 
of the corresponding ideal gas at the same (T, V, n): 


^depart H, n) = F{T, H, n) - E‘g(r, V, n). (31) 

Ideal gases by definition consist of non-interacting molecules of zero volume. The departure function 
is therefore a measure of the contribution to F(T, V, n) from intermolecular forces plus the non-zero 
volume taken up by the molecules. These contributions can be divided into two categories 


^depart ^depart , ^depart 

^ "^phys assoc ’ 


(32) 


where is the contribution from physical (van der Waals) interactions and non-zero molecular 

volumes, and is the contribution from association (hydrogen bonding). Note that because 

intermolecular forces like ion-dipole interactions and ion-ion interactions are not accounted for in 
the departure function, the CPA EOS cannot model mixtures involving brines, such as H 2 O-CO 2 - 
NaCl mixtures. It can model mixtures with CO 2 and H 2 O only. 


B.2.2 Physical contribution to the Helmholtz energy departure function 

The physical contribution is modeled by the cubic part of the EOS. Unlike the original 

CPA introduced by Kontogeorgis et al., Li and Eiroozabadi use the PR EOS, not the SRK EOS, 
for their cubic part. The expression for Ep^Py* given by the PR EOS is 


^depart 

'^phys 

NRT 


= — ln(l — bp) 


2^/2bRT 


In 


1 -|- (1 -|- ■\/2)bp 

1 + (1 - V2)bp\ ’ 


(33) 


where a = a(T, z) is the temperature-dependent parameter that accounts for the physical interac¬ 
tions, and b = 6(z) is the temperature-independent parameter that accounts for non-zero molecular 
volumes. In mixtures, a and b are computed from the pure component properties ai and bi using 
the standard van der Waals mixing rules: 


aij = (1 - kij)ay'^ay‘^, (34) 

C C 

i=l j=l 
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(35) 


C 

b{z) = '^Zibi, 

i=l 

where kij = kij{T) is the binary interaction coefficient between component i and j. This coefficient 
is fitted to experimental VLE data. Fits to binary water/C 02 mixture data give 


kij = 0.5994(r/Tc,co2) - 0.5088, 

for kij between water and CO 2 in which Tc^c 02 is the CO 2 critical temperature. Correlations for 
kij of non-self-associating components may be found in m- The pure component properties a, and 
bi in (34) and (35) are computed from three parameters if it is a non-self-associating: the reduced 
temperature Trj] the reduced pressure Pr^; and the acentric factor w*. They are defined as 


Tr,i — T/Tc,ii 

Pr,i = P/Pc,i, 


UJi = - log 


10 


pvap 

Pr i. 


-1, 


Tr.i=0.7 


where is the critical temperature of component i, Pcj is its critical pressure, and is its vapor 
pressure. For the acentric factor, the vapor pressure is evaluated at a temperature corresponding 
to Tri = 0.7. For a non-self-associating component i |36], its properties a* and bi are 


ai{T) = ac{Tcj)a{u}i,Trj), 
ac(Tc,i) = 0.45724- 


P ’ 

C,2 




l + mi(l - \lTr, 


nii = 


'0.3796 -h 1.485a;i - 0.1644a;2 -h 0.01667wf, 0.1 < w* < 2.0, 

0.37464 -h 1.54226wi - 0.26992a;2, oji < 0.1, 


bi = 0.0778 


RTc,i 

PcA 


For water (the self-associating component), OHaO is computed from the correlation suggested by 
Mathias et al. j52j : 

«H2 o(T") =00 1 + Cl (^1 — + C2[l — \jTr;R2o) + C 3 (^1 — 

where uq, ci, C 2 , C 3 are constants. Li and Firoozabadi have found that choosing ci = 1.7557, C 2 = 
0.003518, C 3 = —0.2746, and oq = 0.9627 L^-bar/mole^ provides a good match with experimental 
data. The volume parameter of water is 6 H 2 O = 0.01458 L/mole. Equation (33) is commonly 
expressed in terms of the compressibility factor Z 


Z = 


PV 


p 


NET pRT' 


(36) 
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by defining new variables A and B that are analogues to a and b, respectively: 


A{T,P,z) = 


B{T,P,z) = 


aP 

hp 


RT' 


Substituting (36|-(38) into (33) yields 


FI 


depart 


phys 


NRT 


= In 



^ In 

Z P (^1 \/T)B 

[z-Bj 

2 ^/2H 

Z + {l-^/T)B_ 


(37) 

(38) 

(39) 


B.2.3 Association contribution to the Helmholtz energy departure function 


As mentioned in Section B.l[ there are two types of association sites: hydrogen-bond donor sites, 
which we will denote as a sites, and hydrogen-bond acceptor sites, which we will denote as /3 sites. 
Hydrogen bonds form only between sites of opposite types; a (/3) sites on a molecule can associate 
with only the /S (a) sites on another molecule. The general expression for is 


^depart 

assoc 

NRT 


= E 


^ Vij 

2=1 j=a,^ 


in Xij - 


Xij -1 


(40) 


where r]ij is the number of j sites (j = a or /3) on a molecule of component i, and Xij is the fraction 
of component i molecules whose j sites are not associated (i.e., j sites that are still free to form 
hydrogen bonds) [13]. We make the following key assumption: water is the only self-associating 
component in the mixture. Other components may cross-associate with water, but they do not 
self-associate nor do they cross-associate with each other. This means that in addition to being 


inapplicable to electrolytic solutions [see discussion below Equation (32)], the CPA EOS that we 
consider is also limited to mixtures where water is the only self-associating component. It would not 
be applicable to aqueous mixtures of alcohols or amines, for example. Nonetheless, this assumption 
is sufficient for CO 2 -H 2 O mixtures, since CO 2 lacks the hydrogen atoms that are necessary for 
self-association. The site fractions for water (denoted by the subscript ‘w’), are obtained by solving 
equations suggested by Kontogeorgis et al. m- 


1 


Xwa — 


I T PX)i=l ^i'nil3Xil3^wa,i 

_ 1 

1 + p X;f=l ZiPiaXiaKp,! 


(41) 


(42) 


in which 5u,a,ii3 is the association strength between an a site on a water molecule and a (3 site 


on a molecule of component i. The quantity is defined similarly. Equations (41) and (42) 


essentially state that the fraction of non-associated a (/3) sites on water, as expressed by Xwa 
iXwp), decreases with the density of available, meaning non-associated, /3 (a) sites on each compo¬ 
nent i, which is given by pziPipXip {pziPiaXia), and with the association strength of /3 (a) sites on 
component i. The association strengths are defined according to 


^wa,wp ^w0,wa 9^ 


exp 


ksT 


- 1 


(43) 
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( 44 ) 


^wlS,ia — SwP,ia^wl3,wa * 7 ^ (^ 5 ) 

where Swa^ifi = Swa^i/siT) is the parameter of cross-association between an a site on water with a /3 
site on component i, and is defined similarly. The cross-association parameters are expressed 

as second-order polynomials of the reduced temperature Tr^i whose coefficients can be fitted to 
experimental data [2] . The energy and volume parameters of water self-association are denoted by 
e and k, respectively. Their values are 


ks 


= 1738.4 K, 


(46) 


K = 1.8015 X 10 


-3 


L 

mole’ 

where is the Boltzmann constant. These parameters have been fitted to data for pure wa¬ 
ter. Finally, g = g{'q) in (431 represents the contact value of the radial distribution function for 
water. We may denote this function with its argument g to distinguish it from the molar Gibbs 
energy g(T,P,z). Li and Firoozabadi [2] use the expression for g(g) from the Starling-Carnahan 
EOS, which states 


, , 1 - 0.5r? 

sW = (TV^’ 


( 47 ) 


^ = TT = 


bp 

4 


B 


(48) 
but we 

continue to use (47) to be consistent with the work of Li and Firoozabadi. Since water is the only 


Studies have found expressions for g{g) simpler than that in (471 also work well 


self-associating component and other components may cross-associate only with water, as stated 
at the beginning of this section, the site fractions for these components are 


1 


Xia — 


Xip = 


1 T P^w'H'wPXwp^w/S ,2a 

1 

1 T P^w'nwaXwabwa,il3 


(49) 


(50) 


Note that in non-associating fluids, the association strengths are all zero so that from (41) 


and (49)~(50|, all site fractions become zero and the CPA EOS reduces to the PR EOS. In summary. 


the association contribution to the Helmholtz energy departure function is given by (40), where 


the site fractions are obtained by solving (41)-(42) and (49)--(50). The association strengths that 


affect the site fractions are calculated using (43)-(48). Substituting (33) and (40) into (32), we 


obtain the full expression for the Helmholtz energy departure function 


^depart 

NET 


= — ln(l — bp) 


2V2bRT 


In 


1 -|- (1 -|- \/2)bp 
1 + (1 - y/2)bp 


+ 11 l^Xij 

7=1 j=a,P 



(51) 
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Section 


B.4 


will derive the expression for the pressure from (51). This form of jg 

for this purpose because P does not appear expl icitly. However, for the purposes of the fugacity 
coefficient </>*, which we derive in the Section B.3 it is more convenient to express in terms 

of the compressibility factor Z = P/pRT. To do so, we define new expressions for the association 
strengths according to 


/C = 


kP 


WOL - 9^ 

exp 


■S a, 2 ^Id ct, 22 ;/3 

i / w 

^wp^ioL^w 

13,wa 

i ^ w 


- 1 


Note that = iP/RT)5i,,ki for all i,j, k, 1. The site fractions are given by 

Z 

Xwa — 


Xw/a — 


^ T ^iVil3Xil3^wa,ip 

z 

^ T X/i=l ^iViaXia^wf},ia 

z 


Xioi — 


Xifi 


Z T ZujTj^p'X'wPfi,ia 

z 


(52) 

(53) 

(54) 

(55) 

(56) 

(57) 

(58) 

(59) 


Z T Ay,Q,^j^ 

In terms of Z, the association contribution to the Helmholtz energy departure function is given 


by (40), where the site fractions are obtained by solving (56)-rt59|. The association strengths that 


affect the site fractions are calculated using (52)-(55) and ([46|)-(48). Substituting (39) and (40) 


into (32), we obtain the full expression for the Helmholtz energy departure function 

.Z + (1 + ■\/2)P 


^depart / ^ 

NRT ~ [z-B 


2V2B 


In 


Z + {l-V2)B 




2=1 j=a,j3 


(60) 


B.3 Fugacity coefficient, enthalpy, density 

In this section, we derive the expression for the fugacity coefficient (j)i from the Helmholtz energy 


departure function (60). The fugacity coefficient is an important quantity in phase equilibria 


(solubility) calculations, as discussed in Section A.l The fugacity coefficient can also be used to 


compute the enthalpy and the density, as described in Section A.3 We start with the expression 
for the differential of the Helmholtz energy F. For systems where external fields like gravity are 
negligible and there is only pressure-volume work. 


dF = -5dr - PdV + ^ pidm, 


(61) 


2=1 
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from which it follows that 




r,y,n. 


This expression, when combined with the definition of in (31), leads to 


F, n) - F, n) _ 1 /\ 

RT RT dm ) 


( 62 ) 


T,V,n, 


The left-hand side of (62) is similar to the left-hand side of (14), except that the ideal gas in (14) 


is at the same pressure (but not same volume) as the real fluid, while the ideal gas in (|^) is 
at the same volume (but not same pressure) as the real fluid. We can denote the volume of the 
ideal gas in (14) as = NRT/P so that its chemical potential is fif{T,P,n) = //)®(r, W®, n). 
Similarly, the pressure of the ideal gas in (62) is = NRT/V so that its chemical potential is 
/i‘®(r, V, n) = /i)®(T, P'S, n). The two different ideal gas states are related by 


/rf(r,y,n)-/rf(r,F^n) = 


rV 


/yig 


dV 


dF'. 


T,n 


From the equality of mixed derivatives in (61), we readily obtain 


dP\ 

dm) 


T,V,n, 


dm 

dV 


T,n 


which we apply to the integral 


rV 


/yig 


ig' 


M. 

dV 


dF' = - 


T,n 


r (T) 

Vig \ dm j 


Ty',n, 


dF' = PTln I ^ j = -PTln Z. 


Thus, 


m{T, F, n) - F, n) = m{T, P, n) - [/xf (T, P, n) - PTln Z 

so that combining with and ( |62[ ), we obtain 

1 f d \ 


In 4>i = 


RT [ dm ) 


-InZ. 


(63) 


T,V,n, 


Since the departure function can be divided into a physical contribution and an association con¬ 
tribution [see Equation (32)], we can also divide the fugacity coefficient into a physical contribution 
and an association contribution according to 


RT 


1 


^ dm j 


— In Z, 


py,n. 


and 


1 / f) pdepart \ 

1 J assoc I assoc j 

~ RT \ dm ) 


(64) 


Ty,n, 
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It can be shown (see for example, [37]) that the fugacity coefficient of the PR EOS is 


In = -\n{Z - B) + 


Bi / B AZ \ 

l3 [z-B ~ + 2RZ - R2 ) 


2 V 2 B 


2 X)_7 = l ^j^ij 

~A 


^ + (1 + y/2)B 
Z + {1-V2)B 


( 65 ) 


where 


* RT’ 


Aij — 


(XijP 


(1 - kij)ay^ay^P 


J12T2 


( 66 ) 


(67) 


Equation (64) is difficult to evaluate because the site fractions x = ixia,Xii 3 ^X 2 a, ■ ■ ■ ,Xc/3) 


are implicit functions of each other [see (56)-59)], and thus, it is not immediately clear how to 
evaluate their derivatives. We avoid having to evaluate the site-fraction derivatives by employing 
the method developed by Michelsen and Hendriks |5^. Their method involves constructing a 
function Q = Q{T,V,n,x) such that Q is equal to at a stationary point (labeled ‘sp’) 

with respect to the site fractions X- That is, Q = Qsp = at a point where 


dQ 


sp 


dx^ 


V 


= 0 for all Xij- 


( 68 ) 


{T,V,n,Xii) 


The subscript Xij fti® derivative above means that all site fractions besides Xij ^.re held fixed. 
Therefore, 


In^assoc ^ 


/ f) pdepart \ 

1 assoc 1 


f dQsp \ 

\ dm j 

T,V,n, 

V dm ) 


T,V,n, 


Applying the chain rule. 


7 dQsp 
V dm 


T,V,n, 


f dQsp 
V dm 


+ E E 

TRriiiX i=l j=a,l3 


dQsp 


(Ty,n,x„-) 


( dXij 
V dm 


By definition of the stationary point in (68), Equation (69) simplifies to 

'dQsp' 


dm 


T,V,n, 


dQ. 

dm 


sp 


T,V,n„x 


T,V,n, 


(69) 


(70) 


The site fractions x held constant in (70). The result of all this is that the troublesome 
derivatives {dxij/dm)\j’y^- do not need to be evaluated. The task remains to select a good choice 
for the function Q. One such choice is 
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(71) 


C 

Q{T,V,n,x) ='^ni ^ rjij {In Xij - Xij +'^) 

i=l j=a,l3 

1 c c 

~ ^77 ^ V ^ ^ ^ ^ ^ VklXkl^ij.kl- 

i=l j=a,(3 k=l 1=0.,p 

It must first be shown that Qsp = -^assoc’^*/-^^- ^1 ^ stationary point, 


5Qi 


sp 


V ) 

Rearranging, we have 


(T,Rn,x„-) 


/I \ ‘ 71' ' - , 

= - 1 -^ 'nklXkl5ij,kl = 0. 

^ fc=l Z =«,/3 


1 C 

^ ^ ) IJklXkl^ijjkl — !• 

^ Xij 


(72) 


fc=l 1 = 0 , f3 

Substituting this result into (0 and comparing with ( |40[ ), we see (since Zi = rii/N) tha'; 

«.p = E"* E ».i (i“»^ - ^) = 1^. 

i=l j=o,/3 ^ ^ 

Substituting into ( [70| ), we have 


dQ, 


sp 


dn^ 


T,V,Ilm,X 


= (InXmi - Xmj + 1) - a,pVmjXmj ^ ^ r]klXklSmj,kl 


J=0‘,l 


2R 


E VijXijY^n-k ^kiXki 


k=l 1=0,0 
f d5ij^i^i 


2=1 j=a,0 


k=l 1=0,0 


V drin 


T,V,nm.,X 


Using (721, this simplifies to 


dQsp \ 
dTlfn ) 


T,V,nm,X 


— E/ 

j = Oi,P 


^ c c 

“ w E E E E ^kixki 

2=1 j=o,0 fc=l 1=0,0 


( dSjj^ki \ 
V y 


T,y,n, 


:X 


From (431, we see that the only part of 6ij^ki that depends on the mole numbers is the contact value 
of the radial distribution function g{r]). Therefore, using (48) yields 


/ ddij^ki 
V dn-m 


T,V,n.m,X 


ij,kl 

/ (9ry \ 


9 dry 

\dnmJ 

Ty,nm,x 


dg _ 

2.5 - ry 


dry 

(1 -ry)4 
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^ij,kl dg bm 

g dry 4U 


bij^kl ^g Bjji 1 


where 





























Substituting into the above, we have 


f dQ: 


sp 


V dur. 


— ^ ^ Vmj 111 Xmj 
T,V,nm,X j=a,fS 


Bm dg 

8gZ drj ^ 

^ ' 1=1 j=a,/3 


E I y E E ^kiXkiS: 


ij,kl 


k=l l=a,P 


Substituting (721 gives 


In^assoc^ ^ Vij'^^Xij + E VkjiXkj-l)- 

ogZ ary 


3 = 0 -, 13 


Combining (65) and (73), we obtain the full expression for the fugacity coefficient 


(73) 


ln(/)j =-ln(Z - i?) + E ^ 


TlZ 


B \Z-B Z‘^ + 2BZ-B‘^ 


A 

2 V 2 B 


2 X)j=i Bi 


A 


B 


In 


■^ + (1 + \/2)B 


Z + {l-V2)B 


Bi dg 


+ E ’!«lnx«+ ^ 

j=a,g ^ ' k=l j=afi 


E^*: E gkjiXkj-^)- 


(74) 


For completeness, we also present the formula for In cjii in terms of the molar density p, which can 


be obtained from (74) using the definitions (36|-(38) and (66|-(67| 


\Ti(f)i = — In (1 — bp) —lnZ+-d^ 


ap 


b \l — bp i7T 1 + 26p — (6p)2 


2y/2bRT \ 


/ 2 X)j=l 


-^r 


1 + (1 + '\/ 2 )bp 

iTTr^E^ 


hipdg 


+ E + E gkj{xkj-i) 


3 = 0,0 




(75) 


0 


B.4 Pressure and compressibility factor 

The main purpose of this section is to derive the expression for the pressure from the Helmholtz 


energy departure function. From (61), the pressure is 


\dV 


T,n 


SO that if = NRT/V is the pressure of the ideal gas at (T, V, n), 

/ a pdepart \ 

p _ pig = _ - 

V Jrn 
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As mentioned in Section B.2.3 


these derivatives are difficult to evaluate because i^depart | 0 Qj 


depends explicitly on pressure, and so the right-hand side will contain explicit volume terms as 
well as terms which depend on the volume implicitly through the pressure. We can overcome this 
problem by noting that since all of the mole numbers in these derivatives are held constant, the 
total number of moles N is also constant. The condition that all mole numbers n be held constant 
is equivalent to the condition that all mole fractions z be held constant. (The converse is not true, 
however.) We can therefore write the pressure in terms of the molar density p = N/V as 




N \ dp J 


T,z 


The pressure and the compressibility factor Z = P/pRT are related to the Helmholtz energy 
departure function according to 


N \ dp I 


T,z 


Z - 1 


1 f d \ 

NRT ap ) 


(76) 


(77) 


T,z 


These expressions are useful because we can use (51) for i^depart^ which is the form of in 

which P does not appear explicitly. 

The pressure is derived from _pdepart ^gjj^g ( 70 |_ Since Ffs|P“*, we can divide 


the pressure into a physical contribution (from the PR EOS) and an association contribution, just 
like we did for the fugacity coefficient in the previous section. The two contributions to P are 


pphys _ pig _|_ 


I / <hT ) 

N \ dp j 


— pRT P 


T,z 


I 

N \ dp j 


N y dp j 


T,z 


The derivative of Fpj^y^* is 


r,z 


(78) 


( OpdepartX 
phys I 

9 . ) 


T,n 


= NRT-^ (-ln(l - bp) -- 

dp \ ^ 2y/2bRT 


In 


1 -|- (1 -|- \/2)6p 

1 + (1 - V2)bp 


T,z 


= NRT 

= NRT 


(1 -|- \/ 2)6 


(l-\/ 2)6 


1 - bp 2V2bRT V1 -h (1 -h y/2)bp 1 -h (1 - \^)V, 
5 1 a 


1 — bp RT 1 + 2bp — (bp)"^ 


so that pP*^y® is 
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pphys 


pRT + p^RT 


b 

1 — bp 


1 a 

RT l + 2bp- {bp)‘^ 


pRT ap^ 

1 — bp 1 + 26/9 — (6/3)2 ' 


( 79 ) 


For the derivative in (781, we encounter the same issue as was encountered with the fugac- 
ity coefficient. Namely, we have to evaluate derivatives of the site fractions. We use the same 
method from Michelsen and Hendriks [53] to overcome this problem. The only difference is that 


Q depends on density instead of volume and on mole fractions instead of mole numbers. That is, 
Q = Q{T, p, z, x). Let us choose Q to be 


Q{T,p,z,x) H (InXii - + 1 ) 

2=1 j=a,l3 

c c 

VijXijJ2^k VkiXkiSij,ki- 
2=1 j=a^l3 k=l 1=C(,I3 

It must first be shown that Qsp = F^^%^^/NRT. At a stationary point. 


(80) 


dQ 


sp 


V ^Xij 

Rearranging, we have 


(Tp,z,Xii) 


= ZiPij I — - 1 j - ZiPij \p^Zk ^ PkiXkAjM 
/ \ k=l l=a,l3 


= 0 . 


1 

p 'y Zk 'y ( PklXkl^ijjkl — 

k=l l=a,P 


Substituting this result into (801 and comparing with (40), we see that 

1=1 j=a,P 


(81) 


Taking the derivative of (80) and using (811, we have 



Equation (43) shows that the only part of bij^ki that depends on the density is g. 
using (48) yields 


Therefore, 


/ dSij^ki \ 

) 


T,z,,x 


^ij,kl ^9 
9 dr/ 



^ij,kl d^ ^ 
g dr/ 4 
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Substituting into the above, we have 


1 / ^ ^depart \ 

^ I assoc I 


NET [ dp J 


T,z 


dQ, 


sp 


dp 


T,z,x 


— 2 ^ ^ ''lijixij 1 ) ^ '^dXij ^ VklXkl^ijM 




/c=l l=a,^ 


Substituting (48) and (81) gives 


1 (dF^\ 

NRT \ dp ) 


T,z 


dQsp 

dp 


T,z,x 


^ ^ ' 1=1 j=a,0 


Multiplying both sides by p^RT, we get the association contribution to the pressure 

r? f d pdepart \ 
oassoc r I ^ assoc \ 


N \ dp J 


T,z 


pRT ( n d^ \ 

i=l j=a,p 


We combine (79) and (82) to write the complete expression for the pressure 

pdg 


P = 


pRT 


ap 


+ 


pRT 


I-bp l + 2bp- (6p)2 ' 2 

Dividing both sides by pRT, we obtain the compressibility factor 

1 1 ap 


^ ' 1=1 3=a,P 


z = 


I-bp RTl + 2bp- {bpY ^2 5 di?) ^ !)• 


Using (37|-(38|, we can also express the above as 


Z = 


z 


Z -B Z2 + 2BZ - ^2 




(82) 


(83) 


(84) 


In the absence of associating components, only the first two terms on the right-hand side of (84) are 


retained. These terms come from the PR EOS. We obtain a cubic polynomial in Z if we rearrange 
the subsequent equation to have all terms on one side. This result shows why the PR EOS is 
said to be a cubic equation of state. If the association terms are non-zero, solving for Z becomes 


more difficult because we do not know a priori the number of roots in (84). As a result, iterative 


root-finding procedures (usually a combination of the bisection method with Newton’s method) 
must be used over a large search interval such as Z € [B 6, lOOOR], where 5 is a small number. 
The lower limit in this interval is only slightly larger than B because B = B(T, P, z) represents the 
closest packing (most compressed state achievable) for a fluid at the conditions (T,P,z). 
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B.5 Simplifications employed in the Li and Firoozabadi model 


In addition to the limitation that water is the only self-associating component in the mixture, which 
we have already discussed in Section B.2.3, Li and Firoozabadi [2] make two more assumptions: 


1. All components have a total of four association sites, with two of each type, so that rjia = 
rjip = 2 for all components i. This four-site model has been shown to work well for water jl2] . 

2. The energetics of association is symmetric between the two types of sites so that Xia = Xip 
for all i. We denote these quantities as Xi = Xia = Xip- 


With these simplifications, (40) becomes 

pdepart 


NRT 


= 4 ^ Zi In Xi - 


2=1 


Xi -1 


As a result, the expressions (51) and (60) for the Helmholtz energy departure function simplify to 


^depart 

NRT 


= — ln(l — bp) 


2V2bRT 


In 


1 -|- (1 -|- \/2)bp 
1 + (1 - y/2)bp 


+ 4^Zi 
2 = 1 




and 


^depart / Z \ 

NRT ^ \Z-b) 


A ^ ^ (1 -|- y/2)B 

2^/2B ^[z + {l- V2)B 


2=1 




respectively. The association strengths, cross-association parameters, and site fractions in (53)-(59) 
are now given by a more tractable set of expressions: 


^W,W — 9^ 


exp 


- 1 


.ksT 

^ 7 “ "^7 

z 


Xw — 


Z + 2 ZiXi^w,i 

z 


Xi = 


Li and Firoozabadi have found that fitting to binary water/CO 2 mixture VLE data gives 


Sw,C02 ~ 0-0529r^QQ2 -|- 0.0404Tr^cO2 ~ 0.0693, 


for the cross-association parameter s^^c 02 between water and CO 2 . 
the fugacity coefficient (jii become 


Equations (74) and (75) for 
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In 4>i 


In {Z -B) + 


Bi / B 

l3 \Z-B ~ Z"^ + 2BZ - B'^, 
^ ^ ^i = l 

2^/2B \ A B J ^ 


+ 41n Xi + 


Bj dg 
2gZ dr/ 


C 


^k{Xk - 1), 

k=l 


Z -\- (\ -\- ^/2)B 
Z + {1-V2)B 


( 85 ) 


\n4>i = 


-ln(l-M-lnZ + |(^ 


ap 


RT l + 2bp- (6p)2 


/ 2 J2j=l 



2y/2bRT \ 

+ 41nxj + zkixk - 1). 

2g dr/ 


1 + (1 + '/2)bp 

1 + (1 - V2)bp 


Equations (83|-(84| simplify to 


n_ pRT _ ap"^ UT (^ .( .-U 

1-bp l + 2bp-{bp)^ ^ V ffdr/j^^*^* 


Z = 


1 


1 


ap 


1-bp i?ri + 26p-(6p)2+^(^ + 5dr/)§^*^^* 


Z = 






Z - B Z2 + 2BZ - B2 




The enthalpy results presented in Sections and are calculated using the equations presented in 
this section. As we demonstrate in those sections, this simplified version of the CPA EOS is sufficient 
for molar enthalpy calculations of mixtures containing CO 2 and H 2 O, as well as enthalpies of pure 
components. It has previously been shown to provide excellent agreement with VLE (solubility) and 
density data [2]. It cannot, however, accurately calculate excess enthalpies of CO 2 -II 2 O mixtures. 
In Section we reason why it fails for excess enthalpy calculations and explain how relaxing the 


two assumptions stated in this section, so that the equations in Sections B.3 and B.4 are followed 
instead, may allow for better agreement with experimental excess enthalpy data. 


C Duan and Sun CO2 activity coefficient model 

In electrolytic solutions, concentrations are usually expressed in terms of molality (m) rather than 


mole fractions, so that (21)-(23| must be recast in terms of molality. In the molality convention, a 


natural choice for the activity coefficient Xi of a solute species i is 


lim 7 j = 1. 

rrij^l 


With 7 j defined in this way, the activity a* and chemical potential pi of i are given by 
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a, = 7*771, = 


fi{T,P,mi) 

/r(r,p,m* = i)' 


Hi{T,P,mi) - fif^{T,P,mi = 1) = -RTlna*. 


( 86 ) 


Here, fl™{T, P, m* = 1) and P, m* = 1) are the fugacity and chemical potential of i, respec¬ 

tively, in an ideal mixture where 777 * = 1. These equations apply to the aqueous phase, which is 
the electrolytic phase. For the C02-rich phase, which is assumed to contain only CO 2 and water, 
Duan and Sun [H] choose the reference chemical potential to be that of the pure component in the 


ideal gas state at a pressure P' = 1 bar. Using (12|-(15|, the chemical potential /7*(r, P, z) of 


component i in the C 02 -rich phase, with pressure expressed in units of bars, is 

/7i(r, P, z) = P' = 1 bar) + RT In ZiP + RT In cj)i{T, P, z). 

For pure CO 2 , this equation simplifies to 


I^C 02 {T, P) = /igo 2 (T, P' = 1 bar) + RT\nP + RT In {T, P). (87) 

Substituting (28) and (29) into the definition (|^ of the CO 2 molar enthalpy of solution Ahsoi, 


A/isoi =Hco2(T,P,mco2 -t 0) - hco^{T,P) 
d 


= -T" 


dT 


f 1^002 (^>-P)r77C02 0)^ 

//^C02(^; P)\ 


_V T j 

P,m \ T J 

P_ 


The subscript m indicates that all species concentrations are held fixed. From (86) and (87), 


d / fJ.c02{T,P,mco2 -t 0) 
dT \ T 


_ ^ = 1)' 
p,m“5ri T 


+ R 


f din 7C02 {T, P, 777C02 0) 
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dT 
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( 1^002 (^, P)\ 

_ ^ 

r/7go^(r,P' = lbar)\ 

+ R 

p 

fdln(^C02{T,P)\ 

dT 

1 T J 

P dT 

V ^ J 

[ dT J 


For a brine solution where the only electrolyte is NaCl, Duan and Sun’s CO 2 activity coefficient 
model is 


ln7C02 = 2Ac02-Na’PNa + 2Ac02-Cl"lCl + Cc02-Na-Cl”7NaPlCl- 


( 88 ) 


The coefficients Aco 2 -Na, Aco 2 -Ch Cc 02 -Na-ci are temperature- and pressure-dependent functions 
that represent intermolecular interactions between CO 2 and the dissolved ions. Equation (88) is 
obtained by taking the derivative, with respect to the molality of CO 2 , of a Pitzer-type excess Gibbs 
energy function [35]. Duan and Sun treat ^^qq^{T,P' = 1 bar) and ACO 2 -CI as being identically 
equal to zero. Using many sources of CO 2 solubility data in pure water and NaCl solutions, they 
fit the functions fi’‘^^{T, P,mco 2 — 1)/-^^) Aco 2 -Na) and Cc 02 -Na-ci to expressions of the form 


37 





























Cl + C 2 T + — + c^T'^ +--+ C6P + C 7 Plnr+ — + +ciirinP, 

T 630 - T ^ '^ r 630-T (630 - T)^ 

where the c’s are constant coefficients. The fugacity coefficient (f)Qo^{T, P) is obtained by iteratively 
solving a nonlinear equation, as described in their paper. In a later study, Duan et al. present a 
series of piecewise expressions (curve fits) for 4 >co 2 can be calculated directly, 

without the need for iterative methods [21] . Although these expressions may increase the efficiency 
of CO 2 solubility computations, they cannot be used for enthalpy computations. We have found 
that for certain conditions, they lead to sharp jumps in the enthalpy. This unphysical behavior 
occurs because the enthalpy is obtained from the fugacity coefficient by taking a temperature 
derivative, so that it should be represented by a smooth function of temperature, rather than a 
series of piecewise expressions. We therefore compute <j)co^{T, P) as described in the earlier study. 
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